Mean-field closure parameters for passive scalar 
turbulence 

J E Snellmani 2, M Rheinhardt^ ^, P J Kapyla^ ^, M J 
Mantere^'^, A Brandenburg^'^ 

1. Department of Physics, Gustaf Hallstromin katu 2a (PO Box 64), FI-00014 
University of Helsinki, Finland 

2. NORDITA, RoslagstuUsbacken 23, SE-10691 Stockholm, Sweden 

3. Department of Astronomy, Stockholm University, SE-10691 Stockholm, Sweden 

E-mail: jan.snellman@helsinki.fi 

Abstract. Direct numerical simulations of isotropically forced homogeneous 
stationary turbulence with an imposed passive scalar concentration gradient are 
compared with an analytical closure model which provides evolution equations for the 
mean passive scalar flux and variance. Triple correlations of fluctuations appearing in 
these equations are described in terms of relaxation terms proportional to the quadratic 
correlations. Three methods are used to extract the relaxation timescales from direct 
numerical simulations. Firstly, we insert the closure ansatz into our equations, assume 
stationarity, and solve for r^. Secondly, we use only the closure ansatz itself and obtain 
Ti from the ratio of quadratic and triple correlations. Thirdly we remove the imposed 
passive scalar gradient and fit an exponential decay law to the solution. 

We vary the Reynolds (Re) and Peclet (Pe) numbers while keeping their ratio at 
unity and the degree of scale separation and find for large Re fair correspondence 
between the different methods. The ratio of the turbulent relaxation time of passive 
scalar flux to the turnover time of turbulent eddies is of the order of three, which 
is in remarkable agreement with earlier work. Finally we make an effort to extract 
the relaxation timescales relevant for the viscous and diffusive effects. We find two 
regimes which are valid for small and large Re, respectively, but the dependence of the 
parameters on scale separation suggests that they are not universal. 



PACS numbers: 47.27.E-, 47.27.tb, 47.40.-x 



1. Introduction 



Fluid flows in astrophysical bodies are most often highly turbulent. Direct modeling 
of such high-Reynolds-number flows is currently impossible. Consequently, greatly 
enhanced diffusivities or modifled diffusion operators are often applied in simulations 
[6]. Such models are challenging in terms of the required computational resources, so 
wide-ranging parameter studies cannot be performed. 

An alternative approach is to separate the large and small scales, and 
derive equations for the former in which correlations of small-scale quantities are 
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parameterized. This is usually referred to as mean-field theory, see e.g. [17, 13, 21, 22]. 
Various schemes have been introduced to close the equations for the correlations of small- 
scale quantities. In astrophysics, the second-order correlation approximation (SOCA), 
and the minimal r approximation (MTA, see e.g. [1, 2]) are widely used. The relevant 
relaxation time in MTA has been determined numerically for passive scalar transport [5] 
as well as for the a effect in mean-field electrodynamics [7, 8] . Another approach, related 
to the MTA, where a relaxation term is invoked to describe the higher order correlations 
has been introduced in [18]. In this 'Ogilvie approach' several nondimensional 
coefficients are invoked to describe physically motivated parameterizations of the higher 
order correlations in terms of relaxation and isotropization terms. This model has been 
applied to different physical setups in order to calibrate the coefficients [9, 16, 14, 10]. 

The validity of the various approximations can in principle be tested by comparing 
with direct simulations in the same parameter regime. In practise this is often not easy 
due to the limited parameter range accessible by the simulations. The starting point for 
such studies has been isotropically forced homogeneous turbulence under the influence 
of rotation [12] and/or shear [24]. In a recent work, the timescales related to diffusion 
and isotropization that appear in the Ogilvie approach have been studied [23]. In the 
present study we extend the work of [23] to the passive scalar case. 



2. Mean-field modelling 

2.1. Ideal case 

Let us consider the transport of a passive scalar under the influence of a turbulent 
fluid motion. For simplicity we assume a homogeneous, incompressible fluid and 
neglect at flrst diffusion and viscous dissipation. Then the governing equations for 
the concentration of the passive scalar, C, and the fluid velocity U read 
r)C 

— = -V-{UC)^-U-VC, (1) 

= - (C7 . V)C7 - -VP + F, V • 17 = 0, (2) 

ot p 

where P is the pressure, p is the constant density and is a forcing function with 
V ■ -F = (with the unit 'force per mass'). Upon introduction of a Reynolds averaging 
procedure, indicated by an ovcrbar, C and U arc decomposed into mean and fluctuating 
parts, C = C + c, U = U + u. The fluctuating flelds, represented by lowercase letters, 
are then governed by 

do — — 

— = - w • VC - L/ • Vc - (w • Vc)', (3) 
ot 

^= -(^.V)L7-(t7-V)w-((n-V)n)'--Vp+/, (4) 

ot p 

where the prime indicates extraction of the fluctuating part, e.g., (uc)' = uc — uc. 
Simplifying further, we stipulate the absence of a mean velocity U and assume that 
the forcing has no mean part, i.e., F = f. In the present case, the goal of mean-fleld 
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modeling consists in deriving a closed equation for the mean concentration C. Prom (1) 
and (3), together with U — we obtain directly 

with the mean density of the passive scalar fiux, J- = cu. So the task of closing 
(5) reduces to representing ^ by the mean concentration C. In the standard mean- 
field approach, (3) is solved for a prescribed fluctuating velocity u, usually under some 
simplifying assumptions which inevitably limit the applicability of the obtained results. 
The solution is employed to express in terms of C. Alternatively, one can abstain from 
deriving such an explicit solution for the fluctuating concentration c and instead strive 
for establishing an evolution equation for which of course again has to be closed in 
the sense that the only variables occurring are the mean quantities C and J- themselves. 
Such an equation is obtained by multiplying (3) with u and (4) with c, summing up 
and averaging, arriving at 

— 1 

— = -itV • iuC) - uV • iuc) -ciiu- V)w) cVp + cf. (6) 

at p 

By virtue of the incompressibihty of the fluid the fluctuating pressure p can be expressed 
by the velocity fluctuations: 

which for an infinitely extended medium and vanishing p at infinity is readily solved by 

( dui I dxj duj I dxi )' {x') 



p 



_P r [dui/dxjduj/dxi) {x') ^2^, 
Att J \x -x'\ 

Now we can conclude that the second, third, and fourth terms on the r.h.s. of (6) are 

quadratic in u and linear in c, hence represent third order correlations. Following [10] 

we introduce here the closure assumption 

1 

-wV • (wc) - c ((w • V)u) cVp = (9) 

P Tq 

with a relaxation time tq. Upon further neglect of the correlation cf, we arrive at 
dTi dC Ti — dC Ti 

— - = -UiUj = -^u^ ' (10) 

Ot OXj Tq OXj Tq 

which governs the evolution of ^. Here, TZij = UiU] stands for the Reynolds stress 
tensor. Since a passive scalar would not act back onto the velocity, TZij can here be 
considered given and the closure is completed. Nevertheless, in the presence of rotation, 
shear, or gravity, (4) contributes quadratic correlations to (6) even in the kinematic 
case; see, e.g., [12, 24, 15]. Equation (10) is similar to the penultimate row of Eq. (53) 
in [10] when replacing temperature perturbation by C and neglecting the buoyancy 
term. Note that for small Tq, that is, fast relaxation, J^j will follow the inhomogeneity 
in (10) almost instantaneously, hence 

Ti ^ -ToTZij (11) 
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and we may interpret TeTZij as a turbulent diffusivity tensor. For a discussion of its 
relationship to traditional mean-field results, see Sect. 2.5. 

To facilitate further comparisons to [10, 4], where an additional evolution equation 
for the mean temperature perturbation variance is derived, we give here an analogous 
equation for = Q, although it is not necessary for completing the closure: 

do — 

^ = -2cw- VC-2c(w- Vc). (12) 

We note in passing that the Q term becomes important in reacting flows [4]. Setting 



-2c(w- Vc) = -Q/t7 (13) 
with another relaxation time Ty, the closed equation reads 
dQ — — Q 

-^^-2J= -VC-- (14) 
dt ^ ^ 

and we have full analogy to the last row of Eq. (53) in [10]. 

Until now we have not constrained the properties of the turbulence, in particular 
we have not required isotropy or homogeneity. For example, inhomogencous turbulence 
could be thought of giving rise to position-dependent relaxation times. However, from 
a strict point of view, the r ansatz (9) is only consistent with an isotropic or uni-axial 
u turbulence where the preferred direction of the latter coincides with the direction of 
JF. Consequently, turbulent properties of u must not change along any other direction. 
The same restrictions should of course hold for the concentration fluctuations c, but this 
is in conflict with the presence of the second preferred direction VC in this turbulence. 
Hence, (9) can be strictly justifled only under very specific circumstances. 

For that reason and for the sake of simplicity, we specify now the mean as horizontal 
average, i.e., as average over all x and y. Consequently, all mean quantities depend on 
z only and only the z component of ^ is relevant. If we further restrict u to have at 
best a z anisotropy, then there is only a single preferred direction, namely that of J- 
and the ansatz (9) is legitimate. The system of mean field equations then simplifies to 

dt~ dz' dt ~ ''dz re' dt~ -^'dz T,- ^ ^ 

where TZzz, and Tj could still depend on z and t. Assuming now further homogeneous 
and statistically stationary fiuctuations u and c and a uniform gradient of C, VC = 
(0, 0, G), a stationary regime of (15) should be given by 

= const. = -TeTZzzG, Q = -2tjTzG = 2rQr^TZz,G^. (16) 

Let us now assume that in a direct numerical simulation (DNS) Eqs. (1) and (2) with 
an appropriately defined forcing / and an imposed uniform G are integrated in time 
until a statistically stationary regime is established. Extracting now all mean quantities 
occurring in (16) from the numerical solution and assuming validity of the model (10) 
and (14), it is obviously possible to determine the crucial relaxation times Tej from such 
runs (method Ml). On the other hand, tqj should of course also obey their defining 
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Figure 1. Method M3 of estimating the timescales tqj from decaying J'z and 
Q normalized by their mean values in the stationary state: Q = {Q)z/{Q)zt, 
T z — z)z/ {J')zt^ where (}j denotes averaging over the variable ^. Dotted vertical 
line: switching off of the imposed mean concentration gradient G. tq = (wrmsfcf)"^ ~ 
dynamical time. Dashed red lines: fit by exponentials in I/tq. 



relations (9) and (13). The third-order correlations ■ (wc), c{{u ■ V)^), cVp and 
c{u ■ Vc) are again accessible in the DNS results and open up an independent path for 
determining the relaxation times (method M2). At the same time, it can also be checked 
to what extent the neglect of cf is justified. 

Another approach to extract Tqj is available from decay experiments, for which, 
after having reached a stationary state in the DNS, the imposed gradient of C is switched 
off. Then, according to (10) and (14), ^ and Q should decay uniformly in space and 
exponentially in time with the increment t^^ and r.^^, respectively, and can be identified 
with the decay rates measured in the DNS (method M3). An exemplification of this 
method is shown in Figure 1. 

The goal of this paper consists in systematically testing the validity of the presented 
closure assumptions for a range of Reynolds and Peclet numbers as well as different levels 
of scale separation. From this we expect hints with respect to the validity of the Garaud 
model of turbulent convection [10]. 
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2.2. Non-ideal effects 

Admitting now diffusion and viscous dissipation, we have to add the term kV'^C with 
the diffusivity k on the r.h.s. of (1) and i/'V^U with the kinematic viscosity u on the 
r.h.s. of (2). Consequently, in the evolution equation (6) for JF, the additional terms 
i/cV^M and KwV^c occur on the r.h.s.. Rewriting their sum as 



lyV'j^, - 2i^^£^ + - iy)u,V^'c (17) 



or more symmetric, as done in [10], as 



i'-\-K^2-^ dui v — K d / dui dc \ 

2 * dxj dxj 2 dxj \ dxj ' dxj J 

does nevertheless not allow a representation entirely by the mean flux. Even in the (very 
particular) case k — v the second terms of (17) and (18) remain. As a skyhook, the 
second and third terms are replaced by the r-ansatz-like expression, —J-jrvt^ although 
they contain second order rather than third order correlations. Analogously, on the 
r.h.s. of (12), diffusion requires a term 



2«;cV2c = reV^g - 2k(Vc)2 (19) 

and the second term is replaced by —Q/t^k- Note that diffusion of T and Q modelled 
in this way is obviously determined by the molecular (or microscopic) diffusivities. 

In astrophysical applications the deviation from ideal conditions is usually small, 
and quantities expressing this smallness are given by the Reynolds and Peclet numbers. 
Re and Pe, which reflect the strength of advection relative to diffusion: 

Re = lirms^/i/, Pe = lirms^/re, (20) 

where ^ is a characteristic scale of the turbulence. We will further make use of the 
Schmidt number Sc = = Pe/Re. 

2.3. Summary of method M2 

So, summarizing all the terms used to determine t^vk and T7«;k from Method M2 (ideal 
and non-ideal, see equations (9), (13), (18) and (19)), we have 

1 

— - = UzV ■ (uc) + c((tt • V)Uz) + -c'VzP + 

V — K, 



+ K)Vc-Vu, —{cV^u.-u.V^c), (21) 



Q 



= 2 (c{u ■ Vc) + k(Vc)2^ . (22) 
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2.4- Scaling of the relaxation times 

For the relaxation times r6,7,yK,ftK some reasonable scaling assumptions are in order 
and we follow essentially the choices of [10]: Tej, belonging to third-order correlation 
terms, are expressed as (CejWrms^i)"^) a-nd t^k,kk, belonging to diffusive second-order 
correlation terms are written as 

{C{u + K)kl/2)-^ and {CKkf)-\ (23) 

respectively. The first of these expressions seems to be appropriate only for — <S 
v + K, hence in general the scaling ansatz should read instead 

((a+«^ + a-«(i^ - K))kiy! (24) 

Here ki — 27r/L is the smallest wavenumber consistent with the box size, L. The 
crucial question is now: are the constants Cqj^i^^^^^^^ universal, at least for a given type of 
turbulence, and in particular, arc they independent of the dimensionless numbers of the 
problem, i.e.. Re and Pe, and of the degree of scale separation? A preliminary answer 
to this question was given in [5], where the timescalcs were found to show a slightly 
increasing trend with increasing scale separation (see their Fig. 4). 

Methods Ml and M3 for determining the relaxation times described in Sect. 2.1 
have now to be modified in the following way: In (16) we have to replace tq by 
nruK/{T6 + t^k) = Tq^k and T7 by T7Ti^^/{t7 + t,^^) = tj,^^. Both methods then deliver 
only these aggregates and we have to employ the different scalings of the relaxation 
times to figure out the individual constants C*. In contrast, method M2 has merely to 
be extended to include also the additional second-order correlations showing up in (18) 
and (19), that is, to use expressions (21) and (22). 



2.5. Comparison with traditional results 

A standard mean-field approach to (1), employing SOCA, that is, neglecting [u ■ Vc)' 
in (3) yields straightforwardly 

Ti{t) = - y Ui{t)uj{t - r) —{t - r) dr (25) 

from which, under the assumption of good temporal scale separation, 

Ht) = - u,{t)u,{t - r) dr —it) = -^ij-Q^, (26) 

can be concluded. K^j = Ui(t)uj{t — r) dr = TcUi(f)uj(t) can readily be identified 
as turbulent diffusivity tensor. Here the correlation time is just defined by the 
last identity. This clearly resembles the result (10) with rg being identified with the 
correlation time r^, the more so as for the validity of both (26) and (10) the relevant 
time parameter has, in a sense, to be small. 

Relaxing the assumption of good temporal scale separation, that is retaining (25), 
we observe the presence of the so-called memory effect [11], that is, the influence of 
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dC I dxj at earlier times t — r on the mean flux at time t by virtue of a convolution. 
Performing a Fourier transform with respect to time, this convolution turns into a simple 
multiplication and we can write 

= -kij{u:)G{uj) (27) 

with a frequency-dependent turbulent diffusivity tensor kij. This quantity is directly 
accessible to the test- field method as described, e.g., in [11]. Based on numerical 
simulations, and without resorting to SOCA, it has been found that for homogeneous 
isotropic turbulence a satisfactory approximation is accomplished already by 

kij{uj) = 5ij ^° (28) 

1 - IWTk 

where kq is the turbulent diffusivity for stationary fields and is independent of cj. A 
slightly better fit is provided by 

= 5,,{1 + K)k,-^^^^-^^^ (29) 

with a constant K. Turning back to the physical space the first approximation (28) is 
equivalent to 

— dT — 

J'+T^^^-KoVC (30) 

or 

f = _^£VC-^. (31) 

dt T„ 

Again, there is striking similarity to (10). Thus by comparing to numerical results 
for r^, a further independent way of checking (9) is provided. The second fit (29) gives 

[11] 

_ dT . 2<92^ . _ . dC 



(1 + K)T + 2r.— + = -(1 + K)koV + r.— J (32) 

indicating the potential importance of higher temporal and mixed temporal/spatial 
derivatives. Note that (31) and (32) are only valid for perfect scale separation in space. 

For the general case of imperfect scale separation both with respect to space and 
time, we refer here to [20] , albeit this work deals with the mean electromotive force of 
MHD rather than with the mean fiux of passive scalar transport. In that work, non- 
locality due to imperfect spatial scale separation shows up in the form of a diffusion 
term r)s^^£ in the evolution equation for £ with a diffusivity rjs occasionally even 
larger than, but of the order of the SOCA estimate of the turbulent diffusivity in the 
high-conductivity limit, r]i = t^u'^.^^^/S. Clearly, this value can be very different from the 
molecular diffusivity. When comparing with the diffusion term for J- identified in (17) 
or (18) where only the microscopic diffusivitics occur we have to state that the Ogilvie 
approach deviates in this respect significantly from what we expect from the traditional 
approach. To reconcile them, possible diffusion terms ~ 'V^J^ and ~ V^Q had to be 
taken into account in the parameterizing ansatzes. 
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2. 6. Significance of method comparisons 

Let us finally discuss what is really "tested" by comparisons of the results of the different 
methods M1-M3. For an incompressible fiuid with the specific conditions of our model 
and under the assumption that no higher than the first order temporal derivative occurs, 
an ansatz for analogous to (16), 



is exhaustive, as 

(i) c and hence is quite generally a linear and homogeneous functional of VC (even 
when the correlation cf is not neglected), and 

(ii) G and hence J^z are spatially constant, both in the statistically steady state and 
during the decay of F^- (That is why the diffusive term oc 'W^J-z is absent.) 

Note, however, that Ki — TZzz as in (16) is an assumption because a contribution 
proportional to VC or dFz/dt can also be provided by the triple correlations, the 
diffusive terms or by cf. 

Since the passive scalar C does not influence the turbulent velocity, the coefficients 
Ki^2 are completely determined by u and hence true constants. Consequently, any 
comparison of the methods Ml, M2, and M3 tests the influence of 

(i) the weak compressibility of the fluid in our simulations, and 

(ii) deviations of the simulated velocity turbulence from homogeneity, isotropy and 
statistical stationarity. 

Note, that the flrst influence can in principle be made arbitrarily small by increasing the 
sound speed Cg in the numerical model, likewise the second by increasing scale separation 
and extending time ranges for averaging. 

A comparison of Ml and M3 tests in addition the justification of 

(i) the neglect of higher temporal derivatives of J-'z, 

(ii) the assumption Ki — TZzz which was employed in calculating tqi,^ — I/K2 from the 
steady state. 

On the other hand, a comparison of Ml and M2 tests again the assumption 
Ki = TZzz and specifically to what extent the neglect of the correlation cf is legitimate. 
This affects tq^^ only, so we expect that here the discrepancies between Ml and M2 are 
more pronounced than in Ty^K. 

With respect to the evolution equation for Q in (15), the same conclusions hold 
true as far as only the steady state and the free decay with G = are taken into 
account. However, considering a general transient as a consequence of switching G 
between two non-zero constants, a richer behavior should appear which perhaps allows 
further reaching conclusions. 



dJ'z 



-KiG - 



(33) 



dt 
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3. Numerical setup 

In order to take benefit of tlie capabilities of tlie Pencil Code :j: we solve instead of 
the incompressible system (1) and (2) the corresponding equations for a compressible, 
but isothermal fluid 

U-VU-VH + f + 2p{V- S{U) + S{U) ■ VH/cl) (34) 
U-VH- clV ■ U, (35) 
V • {CU) + kV^C, (36) 

where we employ the pseudo enthalpy H = Cg In p instead of the density; Cg is the 
constant speed of sound and S{U) the trace-less rate-of-strain tensor Sij = {dUi/dxj + 
dUj/dxi)/2 — dij'V -U/S. Interpretation of the results of such simulations in terms of the 
incompressible model of course requires to keep the Mach number Wrms/ Cs small, typically 
< 0.1. Then it is particularly justified to replace the correlation cVp/p, included in the 
T ansatz, by cV/i. 

The equations are solved by equidistant sixth order finite differences in space and 
an explicit third order time stepping scheme with step size control for stability. The 
computational domain is a cube with dimension (27r)"^ and grid resolutions ranging 
from 32'^ to 256^ according to the requirements raised by the values of the Reynolds 
and Peclet numbers and the forcing wavenumber. Boundary conditions are periodic 
throughout. The fluctuating force / is specifled such that it generates an approximately 
homogeneous, isotropic and statistically stationary fluctuating velocity u. During any 
integration timestep, / is a frozen-in linearly polarized (i.e., non- helical) plane wave 
with a wave- vector which is consistent with the periodic boundary conditions and whose 
modulus is close to a chosen average value kf. The wave amplitude is kept fixed whereas 
the wavcvcctor is randomly changing between time steps and hence / is approximately 
5-correlated in time. For further details, see [3]. 

4. Results 

In general, the parameter space is spanned by the dimensionless numbers s — k{/ki 
(degree of scale separation). Re and Pe, in whose definitions (20) we specify the 
characteristic length i for simplicity by 1/kf. In the following, we will however restrict 
ourselves to Re = Pe, that is, Sc = 1 and leave the more general cases to future work. 
By this, we avoid in particular the complication with the scaling ansatz (23) which 
occurs for v ^ k. For methods Ml and M2, all averaged quantities were derived from 
the simulations by performing in addition to the xy averaging and a temporal averaging 
over an interval in which in particular the correlation (p- was found to be statistically 

X freely available at http://pencil-code.googlecode.com/ 



dU 

~~dt 
dH 

Ih 

dC 
'dt 
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Table 1. Results from methods Ml, M2, and M3 for different values of s = kf/ki. 





V = K, 


Sc = 1. 












s 


Re 




Ml 




M2 


M3 


















1.5 


51.02-552.47 


2.09-2.42 


4.14-4.28 


2.36-2.71 


4.10-6.86 


1.85-2.96 


3.64-5.21 


3 


23.44-694.63 


2.12-2.64 


4.29-5.03 


2.15-2.70 


4.29-14.21 


1.94-2.81 


3.94-4.96 


5 


12.62-415.00 


2.29-2.77 


4.54-5.89 


2.31-2.79 


4.54-41.27 


2.43-2.77 


4.93-5.67 


8 


6.87-256.50 


2.40-2.86 


5.59-9.13 


2.41-2.88 


-21.54-34.81 


2.25-2.85 


5.26-7.70 


10 


5.06-204.34 


2.09-2.87 


6.00-8.06 


2.10-2.93 


-10.43-8.93 


2.24-2.99 


3.44-4.50 



stationary. Statistical errors were estimated by dividing this interval into three equally 
long parts and calculating averages over each of them. These individual averages were 
compared to the average over the whole interval, and the largest deviation was taken 
for the error estimate. 

Wc performed a number of simulations, and extracted tq^k and tt^k using methods 
Ml, M2 and M3. The results are summarized in Table 1, where the several simulations 
are grouped together into different sets by the values used for the scale separation s. 
Within each set the Reynolds numbers were varied by changing u {= n). The timescales 
Tq^k and T-j^f^ listed in Table 1 are also illustrated in Figures 2 and 3, respectively. We 
see that all three methods give quite similar results, with 

T6^^,/to ^2. . .3, and tj^^/tq 7 . . . 10. (37) 

There are some exceptions, however. When s = 1.5, methods M2 and M3 yield 
2 ^ T7kk/to ^ 3, unlike method Ml which stays in range given above. Moreover, 
although method M2 always yields positive values for Tq^^, those for ry^K are sometimes 
negative. The absolute values of the M2 results can also be very high. This is because 
the sum of the correlations used to calculate T7kk may become very small. This problem 
manifests itself mainly in the high Reynolds number runs. 

4-1- Universality of the closure ansatz 

As explained in Sect. 2.2, methods Ml and M3 yield only the quantities Tq,^^ and tyi^^. 
from which the constants Cqj^i,k,kk can be extracted as follows: Recalling the scalings 
(23) we have: 

1 Pi 

= CeUrmsh + C^y^iu + k)^, = CrUrmsh + C^^i^Kkl. (38) 

Multiplying by the viscous time Tyisc — yields 

J- = CesRe + a^, ^ = C^sPe + C,«, (39) 

where a tilde indicates normalization by Tyisc- Hence, when considering I/tq^k and 
1/^7kk as functions of Re (or Pe), the wanted parameters should be obtainable by a 
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Figure 2. Comparison of relaxation timescales rgi/K, normalized by the dynamical 
time To = (wrmsfc/)"^, as functions of the Reynolds/Peclet number from methods MI- 
MS. Different symbols refer to different scale separations s — fcf/fci as indicated in the 
legend. 

linear regression analysis. Figure 4 shows both functions (39) for different values of 
s. Obviously, a linear relation is clearly present both for large and small values of Re, 
but with very different fit parameters for the two ranges. Guided by these functional 
dependencies we hence propose as an alternative for (39) 

^ = a.Re + a. + ^, ^ ^ ^ ^, ^ F(Re), (40) 

and analogously for l/f-ji^f^. This ansatz allows to model linear dependencies on Re both 
for small and large arguments, but with different coefficients: 

F(Re) ^ CgsRe + C^^ for Re ^ oo, (41) 

F(Re) ^ f Ce - -^"j sRe + C,, + for Re ^ (42) 

where the slopes may well be different in sign. The constants in (40) were determined 
by a standard fitting procedure and are given in Table 2. For tqi^k the fit is surprisingly 
good throughout, whereas for f7KK this holds true merely for s = 1.5. For larger scale 
separation the ansatz fails to model the pronounced minima at Re ^ 1 visible in Figure 4 
(right panel). 
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Figure 3. Comparison of relaxation timescales rr^K from methods Ml - M3. Negative 
results are not plotted. For explanations see Figure 2. 







Table 2. Fit parameters of the scaling (40) for the results shown in Fij 


jure 4. 


s 






lO^C^ 












1.5 


0.86 


-121.04 


1.93 


8.02 


0.56 


3.59 


0.36 


-0.20 


3 


1.39 


-45.21 


29.23 


15.5 


0.17 


19.31 


0.56 


15.59 


5 


2.26 


-177.86 


3.47 


4.36 


0.44 


51.31 


0.75 


36.01 


8 


3.06 


-88.07 


12.47 


4.05 


-4.54 


158.69 


0.84 


76.88 


10 


3.50 


-13.25 


25.97 


3.08 


-23.86 


311.38 


0.74 


121.19 



The slopes for high Re, that is, Ce and Cy, show a possible saturation with growing 
scale separation s, but the other coefficients do not. Hence, universality of the ansatz 
(40) with respect to scale separation is questionable. 

4-2. Comparison of methods 

Let us next compare the relaxation timescales rgj^^ and ry^K obtained with the different 
methods in more detail. 
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Figure 4. Relaxation times fgjyK (left) and fy^K (right), normalized to the viscous time 
Tvisc as functions of Re = Pc, for different values of the scale separation s, indicated 
at the curves. Dotted lines with symbols: data from method Ml; solid lines: linear 
fits according to (39), separately for low and high Re. Red dashed: approximation by 
(40) with parameters Cg, Cg, C^^., C^^ from a best fit, see Table 2. 

4.2.1. Methods Ml and MS. Figure 5 shows the ratio of the values of tqi^^ and ry^K 
determined by either of the two methods in dependence on Re and s. The M3 values 
differ from the Ml ones by up to ±30 %. For rg^^K the deviations generally diminish 
with decreasing Re and with growing scale separation, falling beneath 10% for s = S, 10. 
However, neither of these tendencies can be confirmed for tj^k. 

In view of the discussion in Sect. 2.6 it is satisfactory to find improving agreement 
between Ml and M3 with increasing scale separation for which our forced turbulence is 
more and more approaching the desired target of isotropic stationary turbulence. 



4.2.2. Methods Ml and M2. As can be seen in Figure 6, for small scale separations 
s = 1.5,3,5 the values of tqi^^ from both methods coincide fairly. Deviations lie within 
errors. For large s, s = 8, 10, however, we find the deviations grow with falling Re, 
being still within errors around Re = 1. We have to conclude that the neglect of cf has 
its strongest effects for low Re and high s. In contrast, the differences between the ryK^ 
values from methods Ml and M2 are much smaller, reaching a significant magnitude 
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Figure 5. Ratio of the values of tq^k, (left) and ttkk, (right) determined by methods 
Ml and M3. Labels indicate the degree of scale separation s. 
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Figure 6. Ratio of the values of r6i/« (left) and ttkk (right) determined by methods 
Ml and M2. Labels indicate the degree of scale separation s. 



(exceeding errors ) only for the higher s = 8, 10 and Re > 1. A possible reason for this 
is insufficient numerical resolution. 

5. Comments and extensions 



5.1. Alternative scaling 

As an alternative to (38) one might consider 

1 1 



CrUrmskf + C^i^Kkj. (43) 



Then by multiplying with the dynamic time tq = (wrms^f) ^ we arrive at 



= C6 + 



Re 



1 



= C7 + 



Ckk 



(44) 
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Figure 7. Relaxation times tquk (left) and f-jnK (right), normalized on the dynamical 
time To as functions of Re~^ = Pe~^, for different values of the scale separation s, 
indicated at the curves. Dotted lines with symbols: data from method Ml; solid lines: 
linear fits according to (44) , separately for low and high Re. 



where the normahzation is now with respect to tq. Figure 7 shows the same results 
as Figure 4, but now with the altered scaling. Again, a linear fit is viable on each 
curve, but only separately for low and high Re. The corresponding fit parameters can 
be found in Table 3. Unfortunately, an overall fit analogous to (40) does here not work 
satisfactorily. 

5.2. An Ogilvie approach for compressible hydrodynamics? 

One could be tempted to treat the system (36) and (35) in the spirit of the Ogilvie 
approach quite analogously to what was shown in Sect. 2.1. In the ideal case and with 
U — 0, H — we have for the fiuctuating fields 
rhi 

— = -(u-Vuy-Vh + f, (45) 

1^ = - (tt • Vh)' - cy . u, (46) 

and for the mean fiux an analogue to (6), but with the term —cVp/p replaced by 
—c^h. To close the system, an evolution equation for the quantity cV/i seems hence 
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Table 3. 

Alternative fit parameters of the scaling (44) for the results shown in Figure 7. 



s 




small Re 








Re 












c 

^ KK 




C 

^UK 


Cr 


C 

^ KK 


1.5 


0.10 


0.60 


0.064 


0.61 


0.077 


1.25 


0.037 


1.08 


3 


0.091 


0.98 


0.0025 


0.94 


0.14 


0.55 


0.061 


0.97 


5 


0.052 


1.37 


-0.26 


1.23 


0.20 


0.32 


0.088 


0.59 


8 


-2.21 


2.29 


-0.43 


1.62 


0.28 


0.53 


0.11 


0.33 


10 


-5.84 


3.39 


-0.11 


1.88 


0.32 


0.72 


0.12 


0.23 



to be indicated. From (3) and (46) we get 



dcVh 

dt 



- clcW ■ u - cV(u • Vh) - V/iV • (wC) - Vh ■ V(wc), (47) 



d(cVh)^ dh du.- dh dC , ^^«7 , - , , 

^— ^ = - -Uj — - c;c ' - third order terms, 48) 

at oxi ox-i oxi oxi oxiOXi 



where the second order correlation cd'^Uj/dxidxj can only partly be expressed by J-'j. 
The remaining parts could be modelled by a r ansatz as used for the diffusion terms in 
Sect. 2.2, but note that here the "diffusivity" is Cg and we have no argument to consider 
the not properly modelled terms as small. 



6. Conclusions 



The main conclusion to be drawn from the present work is that the time scales used 
to model closure terms in the equations for the mean flux uc and the mean square 
concentration are nearly independent of Re for Re > 10 and also nearly independent 
of the scale separation ratio for k{/ki > 3. Expressed in terms of dynamical times, 
the resulting non-dimensional time scales can be referred to as Strouhal numbers whose 
values are around 3 for the uc closure term and around 7 for the closure term. The 
former value is in good agreement with earlier work using the r approximation [5]. 

Equipped with this knowledge, we may now be better justified in using the closure 
hypotheses discussed here for the quantities uc and c^. On the other hand, as explained 
in the present paper, it is quite clear that these closure hypotheses lack thorough 
justification [19]. One should therefore in future strive to find systematic discrepancies 
from the anticipated scalings. One example that we alluded to in the present paper 
is the inhomogeneous case in which the r approach my break down. Future work in 
that direction seems now to be highly desirable, because in virtually all astrophysical 
applications the turbulence is inhomogeneous or at least anisotropic. 
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